---
title: "SPPQ-23-0075.R2"
author: "Replication Code, Appendix"
format: docx
embed-resources: true
code-background: true
---

```{r preamble, echo = F, warning = F, message = F}
# Clear environment
rm(list=ls())

# Load packages
library(foreign)
library(sandwich)
library(MASS)
library(tidyverse)
library(tidymodels)
library(stargazer); library(knitr)
library(kableExtra)
library(haven)
library(xtable)
library(lme4)
library(corrplot)
library(ggmap); library(sf)
library(ggpubr); 
library(gridExtra)
library(Matching)
library(cspp)
library(urbnmapr)
library(modelsummary)
library(fixest)
options(scipen=7)
options(digits=4)
set.seed(20120630)

# Create clusterMod function for standard errors
clusterMod<-function(model, cluster)
{
  require(multiwayvcov)
  require(lmtest)
  vcovCL<-cluster.vcov(model, cluster)

  coef<-coeftest(model, vcovCL)
  #w<-waldtest(model, vcov = vcovCL, test = "F")
  get_confint<-function(model, vcovCL){
    t<-qt(.975, model$df.residual)
    ct<-coeftest(model, vcovCL)
    cse<-sqrt(diag(vcovCL))
    est<-cbind(ct[,1], cse,ct[,1]-t*ct[,2], ct[,1]+t*ct[,2])
    colnames(est)<-c("Estimate", "Clustered SE","LowerCI","UpperCI")
    return(est)
  }
  ci<-round(get_confint(model, vcovCL),4)
  return(list(coef, ci))
}

# Load in dataset for analysis
data <- read.csv("Women Friendliness Dataset.csv")
uniqcan <- read.csv("unique women candidates.csv")
```

```{r data_wrangling, echo = F, warning = F, message = F}
# Split into Dem, GOP subsamples
datdem <- data %>% dplyr::filter(CandidateDemocrat == 1) 
datgop <- data %>% dplyr::filter(CandidateDemocrat == 0)

# Split by Office Category and Party
govdem <- data %>% dplyr::filter(OfficeCat == "Governor" & CandidateDemocrat == 1)
cabdem <- data %>% dplyr::filter(OfficeCat == "Cabinet" & CandidateDemocrat == 1)
othdem <- data %>% dplyr::filter(OfficeCat == "Other" & CandidateDemocrat == 1)
govgop <- data %>% dplyr::filter(OfficeCat == "Governor" & CandidateDemocrat == 0)
cabgop <- data %>% dplyr::filter(OfficeCat == "Cabinet" & CandidateDemocrat == 0)
othgop <- data %>% dplyr::filter(OfficeCat == "Other" & CandidateDemocrat == 0)

# Split into White, Non-White
whitedat <- data %>% filter(white == 1)
nonwhdat <- data %>% filter(white == 0)

# Create race-and-party subsets
dat.demwh <- data %>% filter(white == 1 & CandidateDemocrat == 1)
dat.demnw <- data %>% filter(white == 0 & CandidateDemocrat == 1)
dat.gopwh <- data %>% filter(white == 1 & CandidateDemocrat == 0)
dat.gopnw <- data %>% filter(white == 0 & CandidateDemocrat == 0)


# To create scatterplots of votepct and wfd by party
dat.wfd <- data %>%
  dplyr::group_by(wfd_ten) %>%
  dplyr::summarise(all.mean = mean(na.omit(votepct)),
            all.sd = sd(na.omit(votepct)),
            all.obs = n())
dat.gop.wfd <- data %>% filter(CandidateDemocrat == 0) %>%
  group_by(wfd_ten) %>%
  dplyr::summarise(gop.mean = mean(na.omit(votepct)),
            gop.sd = sd(na.omit(votepct)),
            gop.obs = n())
dat.dem.wfd <- data %>% filter(CandidateDemocrat == 1) %>%
  group_by(wfd_ten) %>%
  dplyr::summarise(dem.mean = mean(na.omit(votepct)),
            dem.sd = sd(na.omit(votepct)),
            dem.obs = n())

dat.office.dem <- data %>% filter(CandidateDemocrat == 1) %>%
  group_by(wfd_ten, OfficeCat) %>%
  dplyr::summarise(dem.mean = mean(na.omit(votepct)),
                   dem.sd = sd(na.omit(votepct)),
                   dem.obs = n())
dat.office.gop <- data %>% filter(CandidateDemocrat == 0) %>%
  group_by(wfd_ten, OfficeCat) %>%
  dplyr::summarise(gop.mean = mean(na.omit(votepct)),
                   gop.sd = sd(na.omit(votepct)),
                   gop.obs = n())

dat.race.dem <- data %>% filter(CandidateDemocrat == 1) %>%
  group_by(wfd_ten, white) %>%
  dplyr::summarise(dem.mean = mean(na.omit(votepct)),
                   dem.sd = sd(na.omit(votepct)),
                   dem.obs = n()) %>%
  filter(!is.na(white)) %>%
  mutate(white = ifelse(white == 1, "White", "Nonwhite"))

dat.race.gop <- data %>% filter(CandidateDemocrat == 0) %>%
  group_by(wfd_ten, white) %>%
  dplyr::summarise(gop.mean = mean(na.omit(votepct)),
                   gop.sd = sd(na.omit(votepct)),
                   gop.obs = n()) %>%
  filter(!is.na(white)) %>%
  mutate(white = ifelse(white == 1, "White", "Nonwhite"))

# Safe, Competitive, Dangerous state comparisons
dat.dwoc.safe <- data %>%
  filter(avg_votes >= 0.55 & CandidateDemocrat == 1 & white == 0)
dat.dwoc.danger <- data %>%
  filter(avg_votes <= 0.45 & CandidateDemocrat == 1 & white == 0)
dat.dwoc.comp <- data %>%
  filter(avg_votes > 0.45 & avg_votes < 0.55 & CandidateDemocrat == 1 & white == 0)

dat.dww.safe <- data %>%
  filter(avg_votes >= 0.55 & CandidateDemocrat == 1 & white == 1)
dat.dww.danger <- data %>%
  filter(avg_votes <= 0.45 & CandidateDemocrat == 1 & white == 1)
dat.dww.comp <- data %>%
  filter(avg_votes > 0.45 & avg_votes < 0.55 & CandidateDemocrat == 1 & white == 1)

dat.rw.safe <- data %>%
  filter(avg_votes >= 0.55 & CandidateDemocrat == 0)
dat.rw.danger <- data %>%
  filter(avg_votes <= 0.45 & CandidateDemocrat == 0)
dat.rw.comp <- data %>%
  filter(avg_votes > 0.45 & avg_votes < 0.55 & CandidateDemocrat == 0)
```

```{r appendix1, echo = F, warning = F, message = F}
## Summary Table, Identity Variables and Outcome
table.identity <- data %>%
  dplyr::select(votepct, wfd_ten, CandidateDemocrat, white, black,
                latinohispanic, aapi, other, postgrad, governor, cabinet,
                otheroffice)
colnames(table.identity) <- c("Vote Pct", "Full WF Scale", "Cand. Dem", 
                              "Cand. White", "Cand. Black", 
                              "Cand. Latino", "Cand. AAPI", 
                              "Cand. Other Race", "Post-Grad Degree",
                              "Governor", "Cabinet Office", "Other Office")
mean.ti <- sapply(table.identity, mean, na.rm = T)
sd.ti <- sapply(table.identity, sd, na.rm = T)
obs.ti <- apply(table.identity, 2, function(x) sum((!is.na(x))))

tab.comb <- cbind(mean.ti, sd.ti, obs.ti)
colnames(tab.comb) <- c("Mean", "Std Dev", "Obs")

options(width = 60)
kable(tab.comb, caption = "",
      digits = 3)
```

```{r mainmodels, echo = F, warning = F, message = F, include = F}
## Table, Aggregated Main Models with county-clustered SEs
ce1 <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = data,
         se = "cluster", cluster = "fips")

ce2 <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.demwh,
         se = "cluster", cluster = "fips")

ce3 <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.demnw,
         se = "cluster", cluster = "fips")

ce4 <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.gopwh,
         se = "cluster", cluster = "fips")

ce5 <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.gopnw,
         se = "cluster", cluster = "fips")

models_ce1_ce5 <- list("All" = ce1, "White Democrats" = ce2, "NonWhite Democrats" = ce3, 
          "White Republicans" = ce4, "NonWhite Republicans" = ce5)
modelsummary(models_ce1_ce5, output = "markdown",
             fmt_statistic(estimate = 3),
             coef_omit = c(-1, -2, -3, -4, -5, -6, -7, -8, -9, -10, -11, -12,
                           -13, -14, -15, -16, -17),
             vcovBS = "bootstrap", R = 1000, stars = T,
             coef_rename = c("(Intercept)" = "Intercept", "wfd_ten" = "WFI (10)",
                             "wfd_sq" = "WFI Squared", 
                             "CandidateDemocrat" = "Democratic Candidate",
                             "victor" = "Winner", "incumb" = "Incumbent",
                             "postgrad" = "Candidate Post-Grad Degree",
                             "mrp_ideology_mean" = "County Ideology",
                             "pct_dem" = "County Pct Democrat",
                             "female_per" = "County Pct Female",
                             "traditionalistic" = "County Pct Traditionalistic",
                             "moralistic" = "County Pct Moralistic",
                             "cathprop" = "County Pct Catholic",
                             "evanprop" = "County Pct Evangelical",
                             "iwpr_rank" = "County IWPR Rank", 
                             "stateper" = "State Candidate Vote Pct",
                             "TotalVotes" = "County Total Votes (Millions)"),
             gof_map = c("nobs", "r2.adjusted", "F"), align = "lddddd",
             title = "The Effect of Women-Friendliness on U.S. Statewide Executive Woman Candidate County-Level Vote Share, 2010-2019. Robust standard errors reported in parentheses.")
```


```{r modelsummary_2, comment = "", echo = F, message = F, warning = F}
## Match-as-preprocess, race

# Define vars
Y <- data$votepct
Tr <- data$white
X <- data %>% 
  dplyr::select(wfd_AAper, wfd_Hispper, wfd_foriegnb, wfd_medincome, 
                     wfd_college, wfd_bluecollar, wfd_schoolage, wfd_urban, 
                     wfd_nonsouth, wfd_age, CandidateDemocrat, victor, incumb, 
                     postgrad, mrp_ideology_mean, pct_dem, female_per, 
                     traditionalistic, moralistic, cathprop, evanprop, 
                     iwpr_rank, stateper, TotalVotes, Year, governor, cabinet)
ID <- data %>%
  dplyr::select(fips, OfficeCat)
t <- na.omit(cbind(Y, Tr, X, ID))
t_noY <- na.omit(cbind(Tr, X))

# Estimate balance on Tr
tr.predict <- glm(Tr ~ ., data = t_noY)
# summary(tr.predict)

Z <- t[,c(3:5, 7, 9, 11:16, 18, 20:25, 27:28)]

match <- Match(Y=t$Y, Tr=t$Tr, X=t[,c(3:29)], 
               Z=Z, BiasAdjust=T, estimand="ATT", 
               M=1, replace=F, caliper = 2.05)
# summary(match)

postmatch.dat <- as.data.frame(rbind(t[match$index.treated,],
                                     t[match$index.control,])) %>%
  as_tibble()
# tr.postdict <- glm(Tr ~ ., data = postmatch.dat[, -c(1)])
# summary(tr.postdict)

# Create new datasets
postmatch.dat$wfd_ten <- rowSums(postmatch.dat[,c(3:12)], na.rm = T)
postmatch.dat$wfd_sq <- I(postmatch.dat$wfd_ten^2)

match.whd <- postmatch.dat %>% 
  dplyr::filter(Tr == 1 & CandidateDemocrat == 1)
match.nwd <- postmatch.dat %>% 
  dplyr::filter(Tr == 0 & CandidateDemocrat == 1)
match.whr <- postmatch.dat %>% 
  dplyr::filter(Tr == 1 & CandidateDemocrat == 0)
match.nwr <- postmatch.dat %>% 
  dplyr::filter(Tr == 0 & CandidateDemocrat == 0)

## Matched dataset, models
mp1 <- feols(Y ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year), data = postmatch.dat,
         se = "cluster", cluster = "fips")

mp2 <- feols(Y ~ wfd_ten + wfd_sq
         + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = match.whd,
         se = "cluster", cluster = "fips")

mp3 <- feols(Y ~ wfd_ten + wfd_sq
         + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = match.nwd,
         se = "cluster", cluster = "fips")

mp4 <- feols(Y ~ wfd_ten + wfd_sq
         + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = match.whr,
         se = "cluster", cluster = "fips")

mp5 <- feols(Y ~ wfd_ten + wfd_sq
         + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = match.nwr,
         se = "cluster", cluster = "fips")

models_mp1_mp5 <- list("All" = mp1, "White Democrats" = mp2, "NonWhite Democrats" = mp3, 
          "White Republicans" = mp4, "NonWhite Republicans" = mp5)
modelsummary(models_mp1_mp5, output = "markdown",
             fmt_statistic(estimate = 3),
             coef_omit = c(-1, -2, -3, -4, -5, -6, -7, -8, -9, -10, -11, -12,
                           -13, -14, -15, -16, -17),
             vcovBS = "bootstrap", R = 1000, stars = T,
             coef_rename = c("(Intercept)" = "Intercept", "wfd_ten" = "WFI (10)",
                             "wfd_sq" = "WFI Squared", 
                             "CandidateDemocrat" = "Democratic Candidate",
                             "victor" = "Winner", "incumb" = "Incumbent",
                             "postgrad" = "Candidate Post-Grad Degree",
                             "mrp_ideology_mean" = "County Ideology",
                             "pct_dem" = "County Pct Democrat",
                             "female_per" = "County Pct Female",
                             "traditionalistic" = "County Pct Traditionalistic",
                             "moralistic" = "County Pct Moralistic",
                             "cathprop" = "County Pct Catholic",
                             "evanprop" = "County Pct Evangelical",
                             "iwpr_rank" = "County IWPR Rank", 
                             "stateper" = "State Candidate Vote Pct",
                             "TotalVotes" = "County Total Votes (Millions)"),
             gof_map = c("nobs", "r2.adjusted", "F"), align = "lddddd",
             title = "The Effect of Women-Friendliness on U.S. Statewide Executive Woman Candidate County-Level Vote Share, 2010-2019. Robust standard errors reported in parentheses.")
```

```{r competitive, echo = F, warning = F, message = F}
dat.dwoc.safe <- data %>%
  filter(avg_votes >= 0.55 & CandidateDemocrat == 1 & white == 0)
dat.dwoc.danger <- data %>%
  filter(avg_votes <= 0.45 & CandidateDemocrat == 1 & white == 0)
dat.dwoc.comp <- data %>%
  filter(avg_votes > 0.45 & avg_votes < 0.55 & CandidateDemocrat == 1 & white == 0)

# Analyse
tr1 <- feols(votepct ~ wfd_ten + wfd_sq
          + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes, 
         data = dat.dwoc.safe,
         se = "cluster", cluster = "fips")

tr2 <- feols(votepct ~ wfd_ten + wfd_sq
          + victor + postgrad + mrp_ideology_mean
         + pct_dem + female_per
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes, 
         data = dat.dwoc.danger,
         se = "cluster", cluster = "fips")

tr3 <- feols(votepct ~ wfd_ten + wfd_sq
          + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes, 
         data = dat.dwoc.comp,
         se = "cluster", cluster = "fips")

models_tr1_tr3 <- list("Safe Dem" = tr1, "Safe GOP" = tr2, "Competitive" = tr3)
modelsummary(models_tr1_tr3, output = "markdown",
             fmt_statistic(estimate = 3),
             vcovBS = "bootstrap", R = 1000, 
             stars = T,
             coef_rename = c("(Intercept)" = "Intercept", "wfd_ten" = "WFI",
                             "wfd_sq" = "WFI Squared", 
                             "CandidateDemocrat" = "Democratic Candidate",
                             "victor" = "Winner", "incumb" = "Incumbent",
                             "postgrad" = "Candidate Post-Grad Degree",
                             "mrp_ideology_mean" = "County Ideology",
                             "pct_dem" = "County Pct Democrat",
                             "female_per" = "County Pct Female",
                             "traditionalistic" = "County Pct Traditionalistic",
                             "moralistic" = "County Pct Moralistic",
                             "cathprop" = "County Pct Catholic",
                             "evanprop" = "County Pct Evangelical",
                             "iwpr_rank" = "County IWPR Rank", 
                             "stateper" = "State Candidate Vote Pct",
                             "TotalVotes" = "County Total Votes (Millions)"),
             gof_map = c("nobs", "r2.adjusted", "F"), align = "lddd",
             title = "The Effect of Women-Friendliness on U.S. Statewide Executive Woman Candidate County-Level Vote Share, 2010-2019. County-clustered standard errors reported in parentheses.")
```

```{r appendix_eightptwfd, echo = F, warning = F, message = F}
# Testing out alternate metrics
## Table, eight components of WFD with county-clustered SEs
mw81 <- feols(votepct ~ manual_wfd_8 + wfd_8_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = data,
         se = "cluster", cluster = "fips")

mw82 <- feols(votepct ~ manual_wfd_8 + wfd_8_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.demwh,
         se = "cluster", cluster = "fips")

mw83 <- feols(votepct ~ manual_wfd_8 + wfd_8_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.demnw,
         se = "cluster", cluster = "fips")

mw84 <- feols(votepct ~ manual_wfd_8 + wfd_8_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.gopwh,
         se = "cluster", cluster = "fips")

mw85 <- feols(votepct ~ manual_wfd_8 + wfd_8_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(OfficeCat), data = dat.gopnw,
         se = "cluster", cluster = "fips")

models_mw81_mw85 <- list("All" = mw81, "White Democrats" = mw82, "NonWhite Democrats" = mw83, 
          "White Republicans" = mw84, "NonWhite Republicans" = mw85)
modelsummary(models_mw81_mw85, output = "markdown",
             fmt_statistic(estimate = 3),
             coef_omit = c(-1, -2, -3, -4, -5, -6, -7, -8, -9, -10, -11, -12,
                           -13, -14, -15, -16, -17),
             vcovBS = "bootstrap", R = 1000, stars = T,
             coef_rename = c("(Intercept)" = "Intercept", "manual_wfd_8" = "WFI (8)",
                             "wfd_8_sq" = "WFI Squared", 
                             "CandidateDemocrat" = "Democratic Candidate",
                             "victor" = "Winner", "incumb" = "Incumbent",
                             "postgrad" = "Candidate Post-Grad Degree",
                             "mrp_ideology_mean" = "County Ideology",
                             "pct_dem" = "County Pct Democrat",
                             "female_per" = "County Pct Female",
                             "traditionalistic" = "County Pct Traditionalistic",
                             "moralistic" = "County Pct Moralistic",
                             "cathprop" = "County Pct Catholic",
                             "evanprop" = "County Pct Evangelical",
                             "iwpr_rank" = "County IWPR Rank", 
                             "stateper" = "State Candidate Vote Pct",
                             "TotalVotes" = "County Total Votes (Millions)"),
             gof_map = c("nobs", "r2.adjusted", "F"), align = "lddddd",
             title = "The Effect of Women-Friendliness on U.S. Statewide Executive Woman Candidate County-Level Vote Share, 2010-2019. Robust standard errors reported in parentheses.")
```

```{r plot1.office, figures-side, fig.show="hold", out.width="100%", echo=FALSE, warning=FALSE, fig.cap="Distributions of candidate county-level electoral performance across the range of the women-friendliness index. The line reports mean vote share for all candidates at each level of women-friendliness, while the histogram reports the number of candidates at each level."}
plot1.office.dem.out <- ggplot(dat.office.dem, aes(X = wfd_ten)) +
  geom_col(aes(x = wfd_ten, y = dem.obs), size = 1, color = "gray25", fill = "white") +
  geom_line(aes(x = wfd_ten, y = 125*dem.mean), size = 1.5, color = "gray25", group = 1) +
  scale_y_continuous(sec.axis = sec_axis(~./125, name = "")) +
  ylab("Frequency") + xlab("Democratic Candidates") +
  theme(axis.text.x = element_blank()) +
  facet_wrap(vars(OfficeCat))

plot1.office.gop.out <- ggplot(dat.office.gop, aes(X = wfd_ten)) +
  geom_col(aes(x = wfd_ten, y = gop.obs), size = 1, color = "gray65", fill = "white") +
  geom_line(aes(x = wfd_ten, y = 125*gop.mean), size = 1.5, color = "gray65", group = 1) +
  scale_y_continuous(sec.axis = sec_axis(~./125, name = "")) +
  ylab("Frequency") + xlab("Republican Candidates") +
  theme(axis.text.x = element_blank()) +
  facet_wrap(vars(OfficeCat))

grid.arrange(plot1.office.dem.out, plot1.office.gop.out, ncol = 1, nrow = 2)
```

```{r plot_wfd_officetype, echo=FALSE, warning=FALSE, fig.cap="Estimated marginal effects of WFI on candidate performance, matched datasets. Each subgroup plotted with differing color, linetype, and shape. Simulated 90% confidence intervals are reported in vertical bars around the estimated marginal effect. Black indicates the confidence interval is wholly positive or negative.", fig.show="hold", out.width="100%"}
# Show distribution of the types of offices, broken by party and by office
governors_dem <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year), data = govdem,
         se = "cluster", cluster = "fips")

governors_gop <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year), data = govgop,
         se = "cluster", cluster = "fips")

cabinet_dem <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year), data = cabdem,
         se = "cluster", cluster = "fips")

cabinet_gop <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year), data = cabgop,
         se = "cluster", cluster = "fips")

other_dem <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year), data = othdem,
         se = "cluster", cluster = "fips")

other_gop <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year), data = othgop,
         se = "cluster", cluster = "fips")


wfd_output <- rbind(tidy(governors_dem) %>%
                       dplyr::filter(term == "wfd_ten" | 
                                       term == "wfd_sq") %>%
                       mutate(model = "gov_dem"),
                    tidy(governors_gop) %>%
                       dplyr::filter(term == "wfd_ten" | 
                                       term == "wfd_sq") %>%
                       mutate(model = "gov_gop"),
                    tidy(cabinet_dem) %>%
                       dplyr::filter(term == "wfd_ten" | 
                                       term == "wfd_sq") %>%
                       mutate(model = "cab_dem"),
                    tidy(cabinet_gop) %>%
                       dplyr::filter(term == "wfd_ten" | 
                                       term == "wfd_sq") %>%
                       mutate(model = "cab_gop"),
                    tidy(other_dem) %>%
                       dplyr::filter(term == "wfd_ten" | 
                                       term == "wfd_sq") %>%
                       mutate(model = "oth_dem"),
                    tidy(other_gop) %>%
                       dplyr::filter(term == "wfd_ten" | 
                                       term == "wfd_sq") %>%
                       mutate(model = "oth_gop"))

## Predicted effect models, aggregated and matched models
wfd_uniq <- sort(unique(data$wfd_ten))
wfd_squniq <- sort(unique(data$wfd_sq))

govdem_out <- wfd_output$estimate[1] + (2 * wfd_output$estimate[2] * wfd_squniq)
govgop_out <- wfd_output$estimate[3] + (2 * wfd_output$estimate[4] * wfd_squniq)
cabdem_out <- wfd_output$estimate[5] + (2 * wfd_output$estimate[6] * wfd_squniq)
cabgop_out <- wfd_output$estimate[7] + (2 * wfd_output$estimate[8] * wfd_squniq)
othdem_out <- wfd_output$estimate[9] + (2 * wfd_output$estimate[10] * wfd_squniq)
othgop_out <- wfd_output$estimate[11] + (2 * wfd_output$estimate[12] * wfd_squniq)

wfd_plot <- as.data.frame(cbind(wfd_uniq, govdem_out, govgop_out, cabdem_out,
                                cabgop_out, othdem_out, othgop_out)) %>%
  gather(key = "variable", value = "value", -wfd_uniq) %>%
  mutate(variable = c(rep("Dem_Governor", 11), rep("GOP_Governor", 11),
                      rep("Dem_Cabinet", 11), rep("GOP_Cabinet", 11),
                      rep("Dem_Other", 11), rep("GOP_Other", 11)))

# Create confidence intervals
iter <- 5000
sim_govdem_out <- mvrnorm(n = iter, mu = governors_dem$coefficients, 
                          Sigma = vcov(governors_dem))
sim_govgop_out <- mvrnorm(n = iter, mu = governors_gop$coefficients, 
                          Sigma = vcov(governors_gop))
sim_cabdem_out <- mvrnorm(n = iter, mu = cabinet_dem$coefficients, 
                          Sigma = vcov(cabinet_dem))
sim_cabgop_out <- mvrnorm(n = iter, mu = cabinet_gop$coefficients, 
                          Sigma = vcov(cabinet_gop))
sim_othdem_out <- mvrnorm(n = iter, mu = other_dem$coefficients, 
                          Sigma = vcov(other_dem))
sim_othgop_out <- mvrnorm(n = iter, mu = other_gop$coefficients, 
                          Sigma = vcov(other_gop))

marg.effect.data <- data.frame()
for(i in 1:iter){
  ME.govdem <- sim_govdem_out[i,2] + 2*sim_govdem_out[i,3]*wfd_squniq
  ME.govgop <- sim_govgop_out[i,2] + 2*sim_govgop_out[i,3]*wfd_squniq
  ME.cabdem <- sim_cabdem_out[i,2] + 2*sim_cabdem_out[i,3]*wfd_squniq
  ME.cabgop <- sim_cabgop_out[i,2] + 2*sim_cabgop_out[i,3]*wfd_squniq
  ME.othdem <- sim_othdem_out[i,2] + 2*sim_othdem_out[i,3]*wfd_squniq
  ME.othgop <- sim_othgop_out[i,2] + 2*sim_othgop_out[i,3]*wfd_squniq
  d <- data.frame(wfd_uniq, ME.govdem, ME.govgop, ME.cabdem, ME.cabgop,
                  ME.othdem, ME.othgop)
  marg.effect.data <- bind_rows(marg.effect.data, d)
}

marg.effect.data <- marg.effect.data %>%
  group_by(wfd_uniq) %>%
  summarise(lb.govdem = quantile(ME.govdem, 0.025),
            ub.govdem = quantile(ME.govdem, 0.975),
            lb.govgop = quantile(ME.govgop, 0.025),
            ub.govgop = quantile(ME.govgop, 0.975),
            lb.cabdem = quantile(ME.cabdem, 0.025),
            ub.cabdem = quantile(ME.cabdem, 0.975),
            lb.cabgop = quantile(ME.cabgop, 0.025),
            ub.cabgop = quantile(ME.cabgop, 0.975),
            lb.othdem = quantile(ME.othdem, 0.025),
            ub.othdem = quantile(ME.othdem, 0.975),
            lb.othgop = quantile(ME.othgop, 0.025),
            ub.othgop = quantile(ME.othgop, 0.975)
            )

wfd_plot$LB <- as.numeric(paste(c(marg.effect.data$lb.govdem, marg.effect.data$lb.govgop,
                                  marg.effect.data$lb.cabdem, marg.effect.data$lb.cabgop,
                                  marg.effect.data$lb.othdem, marg.effect.data$lb.othgop)))
wfd_plot$UB <- as.numeric(paste(c(marg.effect.data$ub.govdem, marg.effect.data$ub.govgop,
                                  marg.effect.data$ub.cabdem, marg.effect.data$ub.cabgop,
                                  marg.effect.data$ub.othdem, marg.effect.data$ub.othgop)))
wfd_plot$party <- c(rep(c(rep("Dem", 11), rep("GOP", 11)), 3))

wfd_plot <- wfd_plot %>%
  mutate(sig = ifelse(sign(UB) == sign(LB), "significant", "not significant"))

# and plot
pd <- position_dodge(0.75)
override.shape = c(16, 16, 18, 18, 20, 20)
override.linetype = c(1, 3, 1, 3, 1, 3)

plot.wfd.office <- ggplot(wfd_plot, aes(x = wfd_uniq, y = value, 
                                    color = sig,
                                    group = variable,
                                    shape = party,
                                    linetype = variable)) +
  geom_errorbar(aes(ymin = LB, ymax = UB), width = 0.25, position = pd) +
  geom_point(position = pd) +
  geom_hline(yintercept = 0) +
  labs(color = "", linetype = "", shape = "") +
  xlab("WFI Values") + ylab("Marginal Effect") + theme_bw() +
  theme(axis.text.x = element_blank(), legend.position = "none") +
  # guides(color = guide_legend(override.aes = list(shape = override.shape, 
  #                                                 linetype = override.linetype))) +
  scale_shape(guide = "none") + 
  scale_linetype(guide = "none") +
  scale_color_manual(values = c("significant" = "gray10",
                                "not significant" = "gray60")) +
  facet_wrap(facets = "variable", ncol = 3, nrow = 2)

plot.wfd.office
```

```{r reg_output_1, comment = '', message=FALSE, echo=FALSE}
models_governors <- list("Democratic Governors" = governors_dem, 
                         "Republican Governors" = governors_gop)
modelsummary(models_governors, output = "markdown",
             fmt_statistic(estimate = 3),
             vcovBS = "bootstrap", R = 1000, stars = T,
             coef_omit = c(-1, -2, -3, -4, -5, -6, -7, -8, -9, -10, -11, -12,
                           -13, -14, -15, -16),
             coef_rename = c("(Intercept)" = "Intercept", "wfd_ten" = "WFI",
                             "wfd_sq" = "WFI Squared",
                             "victor" = "Winner", "incumb" = "Incumbent",
                             "postgrad" = "Candidate Post-Grad Degree",
                             "mrp_ideology_mean" = "County Ideology",
                             "pct_dem" = "County Pct Democrat",
                             "female_per" = "County Pct Female",
                             "traditionalistic" = "County Pct Traditionalistic",
                             "moralistic" = "County Pct Moralistic",
                             "cathprop" = "County Pct Catholic",
                             "evanprop" = "County Pct Evangelical",
                             "iwpr_rank" = "County IWPR Rank", 
                             "stateper" = "State Candidate Vote Pct",
                             "TotalVotes" = "County Total Votes (Millions)"),
             gof_map = c("nobs", "r2.adjusted", "F"), align = "ldd",
             title = "The Effect of Women-Friendliness on U.S. Statewide Executive Woman Candidate County-Level Vote Share, 2010-2019. Robust standard errors reported in parentheses.")
```

```{r reg_output_2, comment = '', message=FALSE, echo=FALSE}
models_cabinet <- list("Democratic Cabinet" = cabinet_dem, 
                         "Republican Cabinet" = cabinet_gop)
modelsummary(models_cabinet, output = "markdown",
             fmt_statistic(estimate = 3),
             vcovBS = "bootstrap", R = 1000, stars = T,
             coef_omit = c(-1, -2, -3, -4, -5, -6, -7, -8, -9, -10, -11, -12,
                           -13, -14, -15, -16),
             coef_rename = c("(Intercept)" = "Intercept", "wfd_ten" = "WFI",
                             "wfd_sq" = "WFI Squared",
                             "victor" = "Winner", "incumb" = "Incumbent",
                             "postgrad" = "Candidate Post-Grad Degree",
                             "mrp_ideology_mean" = "County Ideology",
                             "pct_dem" = "County Pct Democrat",
                             "female_per" = "County Pct Female",
                             "traditionalistic" = "County Pct Traditionalistic",
                             "moralistic" = "County Pct Moralistic",
                             "cathprop" = "County Pct Catholic",
                             "evanprop" = "County Pct Evangelical",
                             "iwpr_rank" = "County IWPR Rank", 
                             "stateper" = "State Candidate Vote Pct",
                             "TotalVotes" = "County Total Votes (Millions)"),
             gof_map = c("nobs", "r2.adjusted", "F"), align = "ldd",
             title = "The Effect of Women-Friendliness on U.S. Statewide Executive Woman Candidate County-Level Vote Share, 2010-2019. Robust standard errors reported in parentheses.")
```

```{r reg_output_3, comment = '', message=FALSE, echo=FALSE}
models_other <- list("Democratic Cabinet" = other_dem, 
                         "Republican Cabinet" = other_gop)
modelsummary(models_other, output = "markdown",
             fmt_statistic(estimate = 3),
             vcovBS = "bootstrap", R = 1000, stars = T,
             coef_omit = c(-1, -2, -3, -4, -5, -6, -7, -8, -9, -10, -11, -12,
                           -13, -14, -15, -16),
             coef_rename = c("(Intercept)" = "Intercept", "wfd_ten" = "WFI",
                             "wfd_sq" = "WFI Squared",
                             "victor" = "Winner", "incumb" = "Incumbent",
                             "postgrad" = "Candidate Post-Grad Degree",
                             "mrp_ideology_mean" = "County Ideology",
                             "pct_dem" = "County Pct Democrat",
                             "female_per" = "County Pct Female",
                             "traditionalistic" = "County Pct Traditionalistic",
                             "moralistic" = "County Pct Moralistic",
                             "cathprop" = "County Pct Catholic",
                             "evanprop" = "County Pct Evangelical",
                             "iwpr_rank" = "County IWPR Rank", 
                             "stateper" = "State Candidate Vote Pct",
                             "TotalVotes" = "County Total Votes (Millions)"),
             gof_map = c("nobs", "r2.adjusted", "F"), align = "ldd",
             title = "The Effect of Women-Friendliness on U.S. Statewide Executive Woman Candidate County-Level Vote Share, 2010-2019. Robust standard errors reported in parentheses.")
```


```{r fixed_random_effects_analysis, echo = F, warning = F, message = F}
# Distribution of county frequency
count_counties <- data %>%
  group_by(fips) %>%
  summarise(count = n()) %>%
  ggplot(., aes(x = count)) +
  geom_histogram(fill = "black", binwidth = 1) +
  labs(x = "Frequency of County in Dataset",
       y = "")

# Fixed Effects
fe1 <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(fips), data = data,
         se = "cluster", cluster = "fips")

fe2 <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(fips), data = dat.demwh,
         se = "cluster", cluster = "fips")

fe3 <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(fips), data = dat.demnw,
         se = "cluster", cluster = "fips")

fe4 <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(fips), data = dat.gopwh,
         se = "cluster", cluster = "fips")

fe5 <- feols(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + as.factor(Year) + as.factor(fips), data = dat.gopnw,
         se = "cluster", cluster = "fips")


# Random effects
library(lme4)
library(nlme)

re1 <- lmer(votepct ~ wfd_ten + wfd_sq
         + CandidateDemocrat + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + (1 | Year) + (1 | fips),
         data = data)

re2 <- lmer(votepct ~ wfd_ten + wfd_sq
         + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + (1 | Year) + (1 | fips),
         data = dat.demwh)

re3 <- lmer(votepct ~ wfd_ten + wfd_sq
         + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + (1 | Year) + (1 | fips),
         data = dat.demnw)

re4 <- lmer(votepct ~ wfd_ten + wfd_sq
         + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + (1 | Year) + (1 | fips),
         data = dat.gopwh)

re5 <- lmer(votepct ~ wfd_ten + wfd_sq
         + victor + incumb + postgrad + mrp_ideology_mean
         + pct_dem + female_per + traditionalistic + moralistic
         + cathprop + evanprop + iwpr_rank + stateper + TotalVotes
         + (1 | Year) + (1 | fips),
         data = dat.gopnw)

```

```{r distribution_counties_frequency, echo = F, warning = F, message = F, fig.cap = "Histogram showing distribution of county frequency in dataset. The average county is in our dataset seven or eight times, while multimodality in the distribution suggests that many counties are only in the dataset on time, four times, or ten times."}
count_counties
```

```{r modelsummary_fe, comment = "", echo = F, message = F, warning = F}
models_fe1_fe5 <- list("All" = fe1, "White Democrats" = fe2, "NonWhite Democrats" = fe3, 
          "White Republicans" = fe4, "NonWhite Republicans" = fe5)
modelsummary(models_fe1_fe5, output = "markdown",
             fmt_statistic(estimate = 3),
             coef_omit = c(-1, -2, -3, -4, -5, -6, -7, -8, -9, -10, -11, -12,
                           -13, -14, -15, -16, -17),
             vcovBS = "bootstrap", R = 1000, stars = T,
             coef_rename = c("(Intercept)" = "Intercept", "wfd_ten" = "WFI (10)",
                             "wfd_sq" = "WFI Squared", 
                             "CandidateDemocrat" = "Democratic Candidate",
                             "victor" = "Winner", "incumb" = "Incumbent",
                             "postgrad" = "Candidate Post-Grad Degree",
                             "mrp_ideology_mean" = "County Ideology",
                             "pct_dem" = "County Pct Democrat",
                             "female_per" = "County Pct Female",
                             "traditionalistic" = "County Pct Traditionalistic",
                             "moralistic" = "County Pct Moralistic",
                             "cathprop" = "County Pct Catholic",
                             "evanprop" = "County Pct Evangelical",
                             "iwpr_rank" = "County IWPR Rank", 
                             "stateper" = "State Candidate Vote Pct",
                             "TotalVotes" = "County Total Votes (Millions)"),
             gof_map = c("nobs", "r2.adjusted", "F"), align = "lddddd",
             title = "The Effect of Women-Friendliness on U.S. Statewide Executive Woman Candidate County-Level Vote Share, 2010-2019. Clustered standard errors reported in parentheses. Coefficients for county- and year- fixed effects omitted.")
```

```{r modelsummary_re, comment = "", echo = F, message = F, warning = F}
models_re1_re5 <- list("All" = re1, "White Democrats" = re2, "NonWhite Democrats" = re3, 
          "White Republicans" = re4, "NonWhite Republicans" = re5)
modelsummary(models_re1_re5, output = "markdown",
             fmt_statistic(estimate = 3),
             coef_omit = c(-1, -2, -3, -4, -5, -6, -7, -8, -9, -10, -11, -12,
                           -13, -14, -15, -16, -17),
             vcovBS = "bootstrap", R = 1000, stars = T,
             coef_rename = c("(Intercept)" = "Intercept", "wfd_ten" = "WFI (10)",
                             "wfd_sq" = "WFI Squared", 
                             "CandidateDemocrat" = "Democratic Candidate",
                             "victor" = "Winner", "incumb" = "Incumbent",
                             "postgrad" = "Candidate Post-Grad Degree",
                             "mrp_ideology_mean" = "County Ideology",
                             "pct_dem" = "County Pct Democrat",
                             "female_per" = "County Pct Female",
                             "traditionalistic" = "County Pct Traditionalistic",
                             "moralistic" = "County Pct Moralistic",
                             "cathprop" = "County Pct Catholic",
                             "evanprop" = "County Pct Evangelical",
                             "iwpr_rank" = "County IWPR Rank", 
                             "stateper" = "State Candidate Vote Pct",
                             "TotalVotes" = "County Total Votes (Millions)"),
             gof_map = c("nobs", "r2.adjusted", "F"), align = "lddddd",
             title = "The Effect of Women-Friendliness on U.S. Statewide Executive Woman Candidate County-Level Vote Share, 2010-2019. Standard errors reported in parentheses. Coefficients for county- and year- random effects omitted.")
```

```{r maps_prep_last, echo = F, warning = F, message = F}
datdem <- data %>% 
  dplyr::filter(CandidateDemocrat == 1) 
datgop <- data %>% 
  dplyr::filter(CandidateDemocrat == 0)

dem_map <- datdem %>%
  group_by(fips) %>%
  summarise(vote_pct = mean(na.omit(votepct)))
gop_map <- datgop %>%
  group_by(fips) %>%
  summarise(vote_pct = mean(na.omit(votepct)))
all_map <- data %>%
  group_by(fips) %>%
  summarise(vote_pct = mean(na.omit(votepct)))
wfd3_map <- data %>%
  dplyr::select(fips, wfd_ten) %>%
  mutate(wfd_three = case_when(wfd_ten >= 0 & wfd_ten <= 2 ~ 0,
                               wfd_ten == 3 ~ 1,
                               wfd_ten >= 4 & wfd_ten <= 10 ~ 0,
                               .default = NA)) %>%
  group_by(fips) %>%
  summarise(wfd_three = mean(wfd_three, na.rm = T)) %>%
  mutate(wfd_three = case_when(wfd_three == 0 ~ "Not Three WFD Measures",
                               wfd_three == 1 ~ "Three WFD Measures",
                               .default = NA))


# Manage shapefile cleaning
shapefile_counties <- get_urbn_map(map = "counties", sf = T)
shapefile_states <- get_urbn_map(map = "states", sf = T)
```


```{r wfd_three_map, echo = F, warning = F, message = F}
# Summary map, WFI == 3 only highlighted
wfdthree_map <- shapefile_counties %>%
  mutate(county_fips = as.numeric(paste(county_fips, sep = ""))) %>%
  left_join(., wfd3_map, by = c("county_fips" = "fips")) %>%
  ggplot() +
  geom_sf(aes(fill = wfd_three, color = NA, size = 0.05)) +
  geom_sf(data = shapefile_states, fill = NA, color = "black", size = 0.25) +
  scale_fill_manual(values = c("gray80", "gray20"), guide = "legend",
                      na.value = "white") +
  scale_color_manual(values = c("black")) +
  theme(panel.background = element_blank(),
        axis.text.x = element_blank(), axis.text.y = element_blank(),
        axis.ticks = element_blank(), text = element_text(size=8),
        legend.position = "bottom",
        legend.key.size = unit(1.25, 'cm'),
        legend.text = element_text(size=10)) +
  labs(title = "Counties with three measures of WFD index", fill = "") +
  ylab("") + xlab("") +
  guides(color = F, size = F, fill = guide_legend(ncol = 3))
wfdthree_map
```

```{r dem_rep_county_map, echo = F, warning = F, message = F}
# county-wide maps
democrats_map <- shapefile_counties %>%
  mutate(county_fips = as.numeric(paste(county_fips, sep = ""))) %>%
  left_join(., dem_map, by = c("county_fips" = "fips")) %>%
  ggplot() +
  geom_sf(aes(fill = vote_pct, color = NA, size = 0.05)) +
  geom_sf(data = shapefile_states, fill = NA, color = "black", size = 0.25) +
  scale_fill_steps2(low = "gray95", mid = "gray45", high = "gray10",
                      na.value = "white", guide = "colourbar",
                      midpoint = 36.4,
                      breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) +
  scale_color_manual(values = c("black")) +
  theme(panel.background = element_blank(),
        axis.text.x = element_blank(), axis.text.y = element_blank(),
        axis.ticks = element_blank(), text = element_text(size=8),
        legend.position = "bottom",
        legend.key.size = unit(1.25, 'cm'),
        legend.text = element_text(size=10)) +
  labs(title = "Democratic Candidates", fill = "") +
  ylab("") + xlab("") +
  guides(color = F, size = F, fill = guide_colorbar(nbin = 100))

republicans_map <- shapefile_counties %>%
  mutate(county_fips = as.numeric(paste(county_fips, sep = ""))) %>%
  left_join(., gop_map, by = c("county_fips" = "fips")) %>%
  ggplot() +
  geom_sf(aes(fill = vote_pct, color = NA, size = 0.05)) +
  geom_sf(data = shapefile_states, fill = NA, color = "black", size = 0.25) +
  scale_fill_steps2(low = "gray95", mid = "gray45", high = "gray10",
                      na.value = "white", guide = "colourbar",
                      midpoint = 36.4,
                      breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) +
  scale_color_manual(values = c("black")) +
  theme(panel.background = element_blank(),
        axis.text.x = element_blank(), axis.text.y = element_blank(),
        axis.ticks = element_blank(), text = element_text(size=8),
        legend.position = "bottom",
        legend.key.size = unit(1.25, 'cm'),
        legend.text = element_text(size=10)) +
  labs(title = "Republican Candidates", fill = "") +
  ylab("") + xlab("") +
  guides(color = F, size = F, fill = guide_colorbar(nbin = 100))

democrats_map
republicans_map
```

